#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "technique_score.h"

#ifndef USER_FRAME_SHIFT
#define USER_FRAME_SHIFT 1
#endif

/* ---------------- 默认配置 ---------------- */
static const TechniqueScoreConfig DEFAULT_TECHNIQUE_CONFIG = {
    .weight_pitch_derivative  = 0.3f,
    .weight_energy_derivative = 0.3f,
    .weight_timbre_derivative = 0.3f,
    .weight_energy_range      = 0.1f,
    .gamma_derivative         = 1.0f,

    .vibrato_freq_min  = 4.0f,   /* Hz */
    .vibrato_freq_max  = 8.0f,   /* Hz */
    .vibrato_bonus_max = 20.0f,  /* 分 */

    .curve_factor = 8.0f
};

/* ---------------- 内部工具函数 ---------------- */
static inline void compute_derivative(const float* input, float* output, int length)
{
    if (!input || !output || length <= 1) return;
    /* 前向差分 */
    output[0] = input[1] - input[0];
    for (int i = 1; i < length - 1; ++i) {
        output[i] = (input[i + 1] - input[i - 1]) * 0.5f;
    }
    /* 后向差分 */
    output[length - 1] = input[length - 1] - input[length - 2];
}

static double calc_pearson_corr(const float* x, const float* y, int n)
{
    if (!x || !y || n < 2) return 0.0;
    double sum_x = 0.0, sum_y = 0.0, sum_xx = 0.0, sum_yy = 0.0, sum_xy = 0.0;
    for (int i = 0; i < n; ++i) {
        double xi = x[i];
        double yi = y[i];
        sum_x  += xi;
        sum_y  += yi;
        sum_xx += xi * xi;
        sum_yy += yi * yi;
        sum_xy += xi * yi;
    }
    double mean_x = sum_x / n;
    double mean_y = sum_y / n;

    double cov   = sum_xy - n * mean_x * mean_y;
    double var_x = sum_xx - n * mean_x * mean_x;
    double var_y = sum_yy - n * mean_y * mean_y;
    double denom = sqrt(var_x * var_y);
    if (denom < 1e-12) return 0.0; /* 避免除零 */
    double r = cov / denom;
    if (r > 1.0) r = 1.0;
    if (r < -1.0) r = -1.0;
    return r;
}

/* ---------------- 颤音检测与评分 ---------------- */
#ifndef VIBRATO_WINDOW
#define VIBRATO_WINDOW 50
#endif
#ifndef VIBRATO_STEP
#define VIBRATO_STEP 25
#endif
#ifndef ENERGY_ACTIVE_FIXED_THRESHOLD
#define ENERGY_ACTIVE_FIXED_THRESHOLD 0.0001f
#endif

typedef struct {
    float vibrato_rate;            /* 频率 Hz */
    float vibrato_extent;          /* 幅度 (MIDI 半音) */
    float vibrato_stability;       /* 峰值强度(0-1) */
    float vibrato_duration_ratio;  /* 持续性 (0-1) */
    int   has_vibrato;             /* 是否存在颤音 */
} VibratoParams;

/* 自相关检测单窗口颤音参数 */
static void compute_vibrato_params(const float* pitch_buf, const float* energy_buf,
                                   int start_idx, float frame_rate, VibratoParams* out)
{
    memset(out, 0, sizeof(VibratoParams));
    /* ---------- 能量检查 ---------- */
    double e_sum = 0.0;
    for (int i = 0; i < VIBRATO_WINDOW; ++i) {
        float e_db = energy_buf[start_idx + i];
        float e_lin = powf(10.0f, e_db / 10.0f);
        e_sum += e_lin;
    }
    float avg_energy = (float)(e_sum / VIBRATO_WINDOW);
    if (avg_energy <= ENERGY_ACTIVE_FIXED_THRESHOLD) {
        out->has_vibrato = 0;
        return;
    }

    /* ---------- 提取并平滑音高 ---------- */
    float buf[VIBRATO_WINDOW];
    int   valid_frames = 0;
    for (int i = 0; i < VIBRATO_WINDOW; ++i) {
        int idx = start_idx + i;
        float sum = 0.0f; int cnt = 0;
        for (int k = -1; k <= 1; ++k) {
            int pos = idx + k;
            if (pos >= 0 && k + i >= 0 && pos < start_idx + VIBRATO_WINDOW) {
                float v = pitch_buf[pos];
                if (v > 0.0f) { sum += v; cnt++; }
            }
        }
        buf[i] = (cnt > 0) ? (sum / cnt) : 0.0f;
        if (buf[i] > 0.0f) valid_frames++;
    }

    out->vibrato_duration_ratio = (float)valid_frames / (float)VIBRATO_WINDOW;
    if (valid_frames < VIBRATO_WINDOW / 4) {
        out->has_vibrato = 0; /* 有效帧太少 */
        return;
    }

    /* 去均值 */
    float mean = 0.0f; int cnt = 0;
    for (int i = 0; i < VIBRATO_WINDOW; ++i) {
        if (buf[i] > 0.0f) { mean += buf[i]; cnt++; }
    }
    mean = (cnt > 0) ? (mean / cnt) : 0.0f;
    float zero[VIBRATO_WINDOW];
    for (int i = 0; i < VIBRATO_WINDOW; ++i) {
        zero[i] = (buf[i] > 0.0f) ? (buf[i] - mean) : 0.0f;
    }

    /* ---------- 自相关 ---------- */
    int lag_min = (int)(frame_rate / 8.0f + 0.5f); /* 8 Hz */
    int lag_max = (int)(frame_rate / 4.0f + 0.5f); /* 4 Hz */
    if (lag_max > VIBRATO_WINDOW - 2) lag_max = VIBRATO_WINDOW - 2;
    if (lag_min < 2) lag_min = 2;

    float best_r = 0.0f; int best_lag = -1;
    for (int lag = lag_min; lag <= lag_max; ++lag) {
        double num = 0.0, d1 = 0.0, d2 = 0.0;
        for (int i = 0; i < VIBRATO_WINDOW - lag; ++i) {
            float x = zero[i];
            float y = zero[i + lag];
            num += x * y;
            d1  += x * x;
            d2  += y * y;
        }
        if (d1 < 1e-9 || d2 < 1e-9) continue;
        float r = (float)(num / sqrt(d1 * d2));
        if (r > best_r) { best_r = r; best_lag = lag; }
    }

    /* ---------- 参数计算 ---------- */
    out->vibrato_stability = (best_r < 0) ? 0.0f : best_r;
    if (best_lag > 0) {
        out->vibrato_rate = frame_rate / (float)best_lag;
    }

    /* 幅度 (半音) */
    float max_dev = 0.0f;
    for (int i = 0; i < VIBRATO_WINDOW; ++i) {
        float dev = fabsf(zero[i]);
        if (dev > max_dev) max_dev = dev;
    }
    out->vibrato_extent = max_dev; /* 单边幅度 */

    /* 判断颤音存在 */
    out->has_vibrato = (out->vibrato_rate >= 4.0f && out->vibrato_rate <= 8.0f &&
                        out->vibrato_extent > 0.5f && out->vibrato_stability > 0.1f);
}

/* 计算整句颤音加分 (0-20) */
static float detect_vibrato_bonus(const ScoreFeatureBuffers* buffers,
                                  int start_idx, int end_idx,
                                  const TechniqueScoreConfig* cfg)
{
    if (!buffers || !cfg || start_idx >= end_idx) return 0.0f;

    const float* pitch_buf  = buffers->pitch_buffer;
    const float* energy_buf = buffers->energy_buffer;

    int total_len = end_idx - start_idx + 1;
    if (total_len < VIBRATO_WINDOW) return 0.0f;

    int window_cnt = (total_len - VIBRATO_WINDOW) / VIBRATO_STEP + 1;
    float frame_rate = 48000.0f / 1024.0f; /* 与特征提取一致: 46.875 Hz */

    VibratoParams best = {0};
    float max_mag = 0.0f;

    for (int w = 0; w < window_cnt; ++w) {
        int win_start = start_idx + w * VIBRATO_STEP;
        VibratoParams params;
        compute_vibrato_params(pitch_buf, energy_buf, win_start, frame_rate, &params);
        if (params.has_vibrato) {
            float mag = params.vibrato_stability * (params.vibrato_extent + 1e-6f);
            if (mag > max_mag) { max_mag = mag; best = params; }
        }
    }

    if (!best.has_vibrato) {
        printf("[Vibrato] no vibrato detected (frames %d-%d)\n", start_idx, end_idx);
        return 0.0f;
    }

    /* ---------- 评分规则 ---------- */
    float score = 0.0f;

    /* 1. 存在 8 分 */
    score += 8.0f;

    /* 2. 频率 (0-3) */
    if (best.vibrato_rate >= 4.0f && best.vibrato_rate <= 6.0f) {
        score += 3.0f * (best.vibrato_rate - 4.0f) / 2.0f;
    } else if (best.vibrato_rate > 6.0f && best.vibrato_rate <= 8.0f) {
        score += 3.0f * (8.0f - best.vibrato_rate) / 2.0f;
    }

    /* 3. 幅度 (0-3) */
    if (best.vibrato_extent >= 0.5f && best.vibrato_extent <= 1.5f) {
        score += 3.0f * (best.vibrato_extent - 0.5f) / 1.0f;
    } else if (best.vibrato_extent > 1.5f && best.vibrato_extent <= 2.0f) {
        score += 3.0f * (2.0f - best.vibrato_extent) / 0.5f;
    }

    /* 4. 稳定性 (0-2) */
    score += 2.0f * best.vibrato_stability; /* best.vibrato_stability 已在 0-1 */

    /* 5. 持续性 (0-4) */
    if (best.vibrato_duration_ratio >= 0.5f) {
        score += 4.0f * (best.vibrato_duration_ratio - 0.5f) / 0.5f;
    }

    if (score > cfg->vibrato_bonus_max) score = cfg->vibrato_bonus_max;
    if (score < 0.0f) score = 0.0f;

    /* ---- 调试打印 ---- */
    printf("[Vibrato] rate=%.2f Hz, extent=%.2f, stability=%.2f, duration=%.2f, bonus=%.2f\n",
           best.vibrato_rate, best.vibrato_extent, best.vibrato_stability,
           best.vibrato_duration_ratio, score);

    return score;
}

/* ---------------- 主函数实现 ---------------- */
float calculate_technique_score(const ScoreFeatureBuffers* buffers,
                                int frame_idx_start,
                                int frame_idx_end,
                                const TechniqueScoreConfig* config)
{
    if (!buffers || frame_idx_start < 0 || frame_idx_end >= buffers->buffer_size ||
        frame_idx_end - USER_FRAME_SHIFT <= frame_idx_start) {
        return 0.0f;
    }

    /* 使用默认配置 */
    TechniqueScoreConfig cfg = (config) ? *config : DEFAULT_TECHNIQUE_CONFIG;

    /* 有效帧区间 (考虑 USER_FRAME_SHIFT) */
    int idx_start = frame_idx_start + 1; /* 导数需要前后各一帧 */
    int idx_end   = frame_idx_end - USER_FRAME_SHIFT - 1;
    if (idx_end <= idx_start) return 0.0f;

    int n = idx_end - idx_start + 1;

    /* 分配临时缓冲 */
    float *der_pitch_user = (float*)malloc(n * sizeof(float));
    float *der_pitch_base = (float*)malloc(n * sizeof(float));
    float *der_energy_user = (float*)malloc(n * sizeof(float));
    float *der_energy_base = (float*)malloc(n * sizeof(float));
    float *der_timbre_user = (float*)malloc(n * sizeof(float));
    float *der_timbre_base = (float*)malloc(n * sizeof(float));

    if (!der_pitch_user || !der_pitch_base || !der_energy_user || !der_energy_base ||
        !der_timbre_user || !der_timbre_base) {
        free(der_pitch_user); free(der_pitch_base);
        free(der_energy_user); free(der_energy_base);
        free(der_timbre_user); free(der_timbre_base);
        return 0.0f;
    }

    /* ------------ 生成原始序列指针 ------------ */
    const float* pitch_base_buf = buffers->pitch_buffer_base;
    const float* pitch_user_buf = buffers->pitch_buffer;
    const float* energy_base_buf = buffers->energy_buffer_base;
    const float* energy_user_buf = buffers->energy_buffer;
    const float* timbre_base_buf = buffers->timbre_buffer_base;
    const float* timbre_user_buf = buffers->timbre_buffer;

    /* 提取并计算导数 */
    for (int i = 0; i < n; ++i) {
        int base_idx = idx_start + i;
        int user_idx = base_idx + USER_FRAME_SHIFT;

        /* 创建窗口: 前 1, 后 1 */
        float p_base_prev = pitch_base_buf[base_idx - 1];
        float p_base_next = pitch_base_buf[base_idx + 1];
        der_pitch_base[i] = 0.5f * (p_base_next - p_base_prev);

        float p_user_prev = pitch_user_buf[user_idx - 1];
        float p_user_next = pitch_user_buf[user_idx + 1];
        der_pitch_user[i] = 0.5f * (p_user_next - p_user_prev);

        float e_base_prev = energy_base_buf[base_idx - 1];
        float e_base_next = energy_base_buf[base_idx + 1];
        der_energy_base[i] = 0.5f * (e_base_next - e_base_prev);

        float e_user_prev = energy_user_buf[user_idx - 1];
        float e_user_next = energy_user_buf[user_idx + 1];
        der_energy_user[i] = 0.5f * (e_user_next - e_user_prev);

        float t_base_prev = timbre_base_buf[base_idx - 1];
        float t_base_next = timbre_base_buf[base_idx + 1];
        der_timbre_base[i] = 0.5f * (t_base_next - t_base_prev);

        float t_user_prev = timbre_user_buf[user_idx - 1];
        float t_user_next = timbre_user_buf[user_idx + 1];
        der_timbre_user[i] = 0.5f * (t_user_next - t_user_prev);
    }

    /* 计算各导数曲线相关性得分 (0-1) */
    double corr_pitch  = calc_pearson_corr(der_pitch_user,  der_pitch_base,  n);
    double corr_energy = calc_pearson_corr(der_energy_user, der_energy_base, n);
    double corr_timbre = calc_pearson_corr(der_timbre_user, der_timbre_base, n);

    float score_pitch_deriv  = (float)((1.0 + corr_pitch)  * 0.5); /* 0-1 */
    float score_energy_deriv = (float)((1.0 + corr_energy) * 0.5);
    float score_timbre_deriv = (float)((1.0 + corr_timbre) * 0.5);

    /* -------- Gamma 压缩 -------- */
    if (cfg.gamma_derivative > 0.0f && cfg.gamma_derivative != 1.0f) {
        score_pitch_deriv  = powf(score_pitch_deriv,  cfg.gamma_derivative);
        score_energy_deriv = powf(score_energy_deriv, cfg.gamma_derivative);
        score_timbre_deriv = powf(score_timbre_deriv, cfg.gamma_derivative);
    }

    /* 能量最大最小差 */
    float user_e_max = -1e9f, user_e_min = 1e9f;
    float base_e_max = -1e9f, base_e_min = 1e9f;
    for (int i = frame_idx_start; i <= frame_idx_end - USER_FRAME_SHIFT; ++i) {
        float eu = energy_user_buf[i + USER_FRAME_SHIFT];
        float eb = energy_base_buf[i];
        if (eu > user_e_max) user_e_max = eu;
        if (eu < user_e_min) user_e_min = eu;
        if (eb > base_e_max) base_e_max = eb;
        if (eb < base_e_min) base_e_min = eb;
    }
    float diff_user = user_e_max - user_e_min;
    float diff_base = base_e_max - base_e_min;
    float score_energy_range = 0.0f;
    if (diff_base > 1e-6f) {
        float rel_err = fabsf(diff_user - diff_base) / diff_base;
        score_energy_range = 1.0f - rel_err;
        if (score_energy_range < 0.0f) score_energy_range = 0.0f;
    } else {
        score_energy_range = 1.0f; /* 基准差极小, 不惩罚 */
    }

    /* 基础得分 (0-1) 加权平均 (删除综合相似度) */
    float base_score =
        cfg.weight_pitch_derivative  * score_pitch_deriv +
        cfg.weight_energy_derivative * score_energy_deriv +
        cfg.weight_timbre_derivative * score_timbre_deriv +
        cfg.weight_energy_range      * score_energy_range;

    if (base_score < 0.0f) base_score = 0.0f;
    if (base_score > 1.0f) base_score = 1.0f;

    /* 组件调试输出 */
    printf("[Technique detail] gamma=%.2f, pitch_deriv=%.3f, energy_deriv=%.3f, timbre_deriv=%.3f, energy_range=%.3f\n",
           cfg.gamma_derivative, score_pitch_deriv, score_energy_deriv, score_timbre_deriv, score_energy_range);

    /* 颤音加分 */
    float vibrato_bonus = detect_vibrato_bonus(buffers,
                                               frame_idx_start + USER_FRAME_SHIFT,
                                               frame_idx_end   + USER_FRAME_SHIFT,
                                               &cfg);

    /* 最终分数 (百分制) */
    float total_percent = base_score * 100.0f + vibrato_bonus;
    if (total_percent > 100.0f) total_percent = 100.0f;

    /* 转换回 0-1 范围 */
    float total_raw = total_percent / 100.0f;

    /* 调试输出 */
    printf("Technique stats: base=%.3f, vibrato_bonus=%.2f, total=%.3f (raw)\n",
           base_score, vibrato_bonus, total_raw);

    free(der_pitch_user); free(der_pitch_base);
    free(der_energy_user); free(der_energy_base);
    free(der_timbre_user); free(der_timbre_base);

    if (total_raw < 0.0f) total_raw = 0.0f;
    if (total_raw > 1.0f) total_raw = 1.0f;
    return total_raw;
} 